Development of a feeding simulation to evaluate how feeding distribution in aquaculture affects individual differences in growth based on the fish schooling behavioral model

In this study, we developed a feeding simulation using the fish schooling behavior model to evaluate growth differences from the feeding spatial distribution. In the proposed simulation, feeding behavior was modeled using the fish schooling model to simulate the amount of feed consumed by each individual. Next, body mass growth was calculated based on the amount of feed consumed. A 3.0-m diameter aquaculture tank was used for the simulation. We used three feeding methods to evaluate how feeding distribution affected growth: Feeding A, B, and C. The feed was distributed in a square pattern with one side length of 1.5, 1.0, or 0.5 m for Feeding A, B, and C groups, respectively. The results revealed that individual differences in body mass resulting from each feeding method differed greatly. The individual difference was largest in the Feeding C group. Here, maximum swimming speed was assumed to be proportional to total length. The feeding area of Feeding C was narrow; therefore, the first individual to arrive in the feeding area dominated the feed. Large individuals accessed the feed more easily than did small individuals. Consequently, the growth of large individuals became more rapid, and the individual differences became large in Feeding C. A rearing test can be conducted in a short time, and the optimal aquaculture operation was easily determined using the proposed simulation method. We concluded that the proposed simulation is useful as a decision-making tool for aquaculture management.


Introduction
The demand for fish has grown dramatically in recent years. Fisheries production from aquaculture has been increasing and exceeded production from fishing vessels in 2015 [1]. Feed is an important aquaculture management component. Less feed leads to less growth of reared fish, whereas overfeeding leads to feed waste. Therefore, a moderate feeding amount is required to minimize economic loss. The FAO reported that feeding cost is the main cost component on aquaculture farms, accounting for 60% of aquaculture production costs [2]. The relationship between the feeding amount and the growth of reared fish has been reported in many studies to determine the optimal amount to feed [3]. Another study reported that the competition for feed is an important factor when rearing fish over a long period [4]. An individual with high competitive ability preferentially accesses food resources; therefore, a highly competitive individual has a high growth rate, and vice versa. Consequently, the competition for feed leads to individual differences in growth. Individual differences in growth increase variability and instability therein. Rearing experiments can be conducted to determine the optimum feeding operation. For example, Handeland et al. [5] performed a rearing experiment targeting Atlantic salmon, Salmo salar, at different water temperatures to determine the water temperature that maximizes feed conversion efficiency, FCE, defined as the ratio of feed intake to body mass. Similar studies have been widely conducted [6]. Other studies have focused on the spatial distribution of feeding. For example, Jørgensen et al. [7] conducted a rearing test in tanks with standing and flowing water. The feed was distributed evenly by water flow in the tanks with flowing water, whereas the feed clumped in the tanks with standing water. The authors concluded that the feed intake variability was small in the tank with flowing water. Therefore, the feeding method, such as the amount and its distribution, affected the final growth of fish.
We hypothesized that the growth difference from feeding in a spatial distribution was determined by the spatial behavior of individual fish in the tank. To understand the growth difference caused by feeding with a particular spatial distribution, the growth history of each individual should be quantitatively evaluated. However, growth is generally evaluated by the mean body mass of reared fish, and growth of an individual fish is difficult to track. Furthermore, it is costly and time-consuming to conduct an actual rearing experiment, and a comparison of various feeding methods is difficult.
Simulation studies have been widely performed on the growth of reared fish to determine feeding efficiency without performing an actual rearing experiment. In some studies, the dynamic energy budget (DEB) theory [8] has been employed to simulate fish growth. For example, Alver et al. [9] applied DEB theory to larval stage cod, Gadus morhua. Føre et al. [10][11][12] conducted a series of simulation studies and developed a simulation model that included the fish behavioral model and DEB theory; they reported the simulated mean growth in body mass, but did not discuss the relationship between individual behavior and individual growth.
In this study, we developed a feeding simulation method based on the fish schooling model to clarify the relationship between individual fish behavior and growth. First, we describe our proposed simulation method and the feeding spatial distribution results. Next, we focus on the spatial distribution of the feed, which reportedly affects the growth of fish. Fish growth was simulated under several feeding spatial distribution conditions using the proposed simulation method. We hypothesized that growth differences according to the feeding spatial distribution were determined by the behavior of individual fish. Finally, we discuss the utility of the proposed simulation method for aquaculture farm management, including the feeding strategy.

Schematic of the proposed simulation model
Rainbow trout, Oncorhynchus mykiss, was targeted in this study. Fig 1 presents a schematic of the proposed simulation model. The proposed simulation was composed of a "feeding simulation" and a "growth simulation" section.   Flow chart of the simulation. For the feeding simulation, the input variables in the growth simulation section were W i n and TL i n . The feeding simulation was conducted for 1,000 s. Finally, S i n was calculated and exported to the growth simulation. The calculated S i n was input into the growth simulation, which calculated W i n and TL i n and returned the values to the feeding simulation. This was repeated until the final day of the experiment (day 90). To confirm the spatial interaction between the fish and feed, a simulation in which randomization was performed every day after the feeding simulation was performed (see "Entire simulation set-up"). the weight of feed consumed by each individual, S i n , was the output. Here, S i n is the weight of feed consumed by the i-th individual on the n-th day. The simulated feed intake weight, S i n , was transferred to the growth simulation section. In the growth simulation section, body mass, W i n , and total length, TL i n , were calculated based on the weight of feed consumed, S i n . The body mass and total length of each individual was transferred to the feeding simulation section. This simulation was repeated until the final day. Unity software 2018.3.12.f1 was used for the simulation and visualization.

Feeding simulation
Schematic of the feeding simulation. Fish behavior, including movement against the feed, was simulated in the feeding simulation section. The weight of the feed consumed was output from the fish behavioral model.
Fish schooling behavior is determined by interference from other individuals. Several models have been developed to evaluate fish school movements from individual behavior, such as the boid model [13], the Sannomiya model [14,15], and the stochastic model [16]. In this study, we used the boid model for the feeding simulation. The boid model was originally proposed to simulate bird flocking behavior [13], and has been applied to fisheries [17].
Here, the motion of an individual fish is expressed using Newton's equation of motion as follows: where W i n is fish body mass and a t , v t , x t , and F t are the acceleration, velocity, position, and swimming force vectors at time, t, respectively. Δt is the time step of the simulation and is set to 0.025 s. F t in Eq (1) is determined based on the boid model [13], in which schooling behavior is determined by an individual's movement affected by neighboring individuals. Fig 3 presents an overview schematic of a modeled individual fish. The individual F t is determined by the other individuals within the field of view. Here, the field of view is assumed to be a simple sphere with a radius of 2TL, except that behind within an angle of 30˚ [17,18]. We assumed that an individual can detect the feed regardless of the field of view. Fig 4 provides an overview of swimming force. Here, fish schooling behavior is expressed by six rules: 1. Separation, F 1 ; 2. Cohesion, F 2 ; 3. Alignment, F 3 ; 4. Avoiding a boundary, F 4 ; 5. Approaching feed, F 5 ; and 6. Random movement, F 6 . As a result, swimming force is calculated using the following equation: where F n is the n-th force vector. The details of each force are explained as follows. Separation. The separation force vector, F 1 , is a repulsive force from the nearest individual, described by the following equation:  where x i is the position vector of the individual being considered, x j is the position vector of the nearest individual, and w 1 is the weight coefficient of F 1 .
Cohesion. The cohesion force vector, F 2 , is defined as the vector that moves toward the center of other individuals in the field of view. As a result, the cohesion force vector is described using the following equation: where x i is the position vector of the individual being considered, and w 2 is the weight coefficient of F 2 . � x is the center of gravity position vector of other individuals in the field of view, calculated using the following equation: where x k is the position vector of the k-th individual in the field of view and m is the number of other individuals in the field of view. Alignment. The alignment force vector, F 3 , is the vector that moves in the same direction as other individuals, and is expressed using the following equation: where v i is the velocity vector of the individual being considered, and w 3 is the weight coefficient of F 3 . � v is the average velocity vector of other individuals in the field of view, calculated using the following equation: where v k is the velocity vector of the k-th individual in the field of view, and m is the number of other individuals in the field of view. Avoiding a boundary. The avoiding boundary vector, F 4 , is the repulsive force from the boundary, including a wall or bottom of the tank and water surface, and is expressed as follows.
where x i is the position vector of the individual being considered, x bound is the position vector of the nearest point of the boundary in the field of view, and w 4 is the weight coefficient of the avoiding force from the boundary. Approaching feed. The approaching feed vector, F 5 , is the attractive force to the feed. The force direction vector is expressed as follows: where x i is the position vector of the individual being considered, x f is the position vector of the feed, and w 5 is the weight coefficient of the feeding behavior force vector. In this simulation, we assumed that the individual can detect the feed regardless of the field of view.
Random movement. The force vector of random movement, F 6 , is considered based on uniformly distributed random numbers as follows: w 6 ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi where w 6 is the weight coefficient of the random movement force, and x 1 , x 2 , and x 3 are random numbers uniformly distributed over the range [-1, 1].

Definition of feed intake during the feeding simulation
In this simulation, contact between the fish model and the feed model was defined as feed intake. The number of feed pellets contacted by each individual, N i f , was counted during the feeding simulation. The contacted feed was deleted immediately during the simulation. As a result, feed intake weight S i n was calculated after the feeding simulation using the equation: where S i n is the feed intake weight of the i-th individual on the n-th day and W f is the weight of a grain of feed. The feed was assumed to be a simple sphere with diameter of 5.0 mm and a volume of 0.065 cm 3 . The density of the feed was the same as water; therefore, weight, W f , was 0.065 g. The sinking speed of the feed was set to 0.025 m/s.

Parameters for the feeding simulation
In the simulation, fish behavior was divided into the "standard mode" and the "feeding mode". Table 1 presents a schematic of the two swimming modes. If the feed was in the tank and the feed intake weight S i n was smaller than the maximum feed intake weight S max , the individual was in "feeding mode". Otherwise, the individual was in "standard mode".
Some studies have investigated feed intake per day; however, we found no study focusing on maximum feed intake weight. Morales et al. [19] reported that the feed intake weight per day was 1.43-2.26% of fish body mass. In this study, we classified excess feeding based on the literature [20,21]; the maximum feed intake, S max was 4% of the body mass (S max ¼ 0:04W i n ). The strength of each force was determined by the weight coefficients, w n described as above. Table 2 lists the w n values. The weight coefficients except for the approaching feed, w 5 , were fixed. However, w 5 was determined by the mode of the individual. Here, w 5 was set to 3.0 for the "feeding mode" and 0.0 for the "standard mode".
Swimming force was calculated using Eq (2), and acceleration, velocity, and position were determined using Eq (1). In this simulation, the maximum swimming velocity magnitude |v| max was defined to consider the limit of fish swimming ability. If the velocity magnitude in Eq (1), |v|, exceeded the maximum velocity magnitude, the swimming velocity vector was corrected as follows: where e v is the unit vector of the swimming velocity vector v. Maximum swimming velocity, |v| max was proportional to total length, TL, as follows: where C v is the coefficient of the maximum velocity magnitude; Table 2 lists the values. The C v value was modified by the swimming modes: C v = 1.5 for the "standard mode" [22] and C v = 7.0 for the "feeding mode" [23].

Growth simulation
The amount of growth of each individual was determined by feed intake, S i n , based on the feeding simulation. We calculated the amount of body mass growth from the feed conversion efficiency, FCE = body mass gain/weight of feed intake, as follows: where W is the body mass, ΔW is the daily body mass gain, and S is the weight of the feed consumed. The subscripts i and n indicate the individual number and the rearing day, respectively. FCE is feed conversion efficiency; here, we assumed FCE to be 1.0 [19]. After the body mass was updated, total length TL i n was calculated from the following allometric equation: where TL i n is the total length of the i-th individual on the n-th day and a and b are coefficients of allometry; here, a was 0.0209 and b was 2.843, respectively [24].

Entire simulation set-up
A 3.0-m diameter aquaculture tank with a depth of 1.5 m was used for the simulation (Fig 5). We simulated the three feeding methods named Feeding A, B, and C to evaluate the impact of the feeding method on growth, particularly the feed distribution (Fig 5). Feed was distributed in a square pattern for all methods, but the one-side length was different. The lengths were 1.5, 1.0, and 0.5 m for Feeding A, B, and C groups, respectively. The feed was homogeneously and randomly distributed using the "Random" function of Unity. In total, 100 individuals were simulated in each feeding method. The feeding simulation was conducted for 1,000 s, which was sufficient time to finish the feeding. The initial position of individuals on day 1 was randomized using the "Random" function of Unity software. After day 2, the initial position was set according to the final position the day before. We assumed that the fish were well dispersed after the feeding simulation. However, as the spatial interaction between the fish and feed is considered important, we then conducted a simulation in which randomization was performed every day.
The growth simulation was initiated after the feeding simulation. In the growth simulation section, the body mass, W i n , and total length, TL i n , were calculated using Eqs (14)- (16). After the growth simulation, the feed intake weight of each individual, S i n , was 0. The entire simulation was continued until day 90. The initial body mass and feeding condition were determined based on a previous growth experiment [19]. Here, the initial body mass of all individuals was set to 38.53 g, and the feed intake weight was 1.86% of body mass per day [19]; therefore, the amount of feed distributed per day was 1.86% of the total body mass of 100 individuals. The feed pellet weight W f was 0.065 g, as mentioned above. The number of pellets given was determined as follows: Amount of feed/W f . The feeding frequency was once per day. The individuals formed a school and swam along the wall of the tank (Fig 6). When feeding started (0-10 s), individuals approached, gathered, and ingested the feed. All of the feed was ingested within 30 s. After finishing feeding (40-50 s), the individuals continued standard swimming, and swam near the wall.

Fish behavior in the tank
The individuals in Feeding C swam along the wall of the tank similar to the Feeding A group (Fig 7). When feeding started (0-10 s), the individuals approached the feed. The individuals gathered in a narrower area compared with that shown in Fig 6. All feed was ingested within 20 s, and the individuals went back to standard swimming.  [19] are shown.

Growth of body mass during the simulation
Body mass increased exponentially in all feeding methods, but the individual difference in body mass increased with the rearing day despite the same initial body mass. In particular, the individual differences in the Feeding C group were very large, and the body masses of several individuals in Feeding C increased dramatically. The mean body mass from the simulation agreed well with the experimental values from Morales et al. [19] for all feeding methods, but the standard deviations were different. The standard deviation of the Feeding A group was smaller than that of the other feeding groups, and the standard deviation increased as time elapsed in the Feeding C group.
Furthermore, as mentioned in the "Materials and Method", after the feeding simulation, we conducted a simulation in which randomization was performed every day (Fig 9). The individual difference shown in Fig 9 is smaller compared with that in Fig 8; some remarkably large individuals are shown in Feeding C in Fig 8, but not in Fig 9. However, the tendency for the individual difference to increase with decreasing feeding area was similar between the two figures. Therefore, the following analysis and discussion are based on the results shown in Fig 8.  Fig 10 presents the body mass distributions on the final day (day 90) of each feeding method. The mean ± standard deviation values of the final body mass were 210.4 ± 15.7 g for Feeding A, 214.5 ± 27.1 g for Feeding B, and 219.4 ± 59.6 g for Feeding C. The skewness of each distribution was 0.180, 0.914, and 1.700 for Feeding A, B, and C, respectively. The means of body mass were slightly but not significantly different among the three feeding methods (Kruskal-Wallis test, p = 0.106). However, a significant difference was detected among the variances of the three distributions (Bartlett test, p < 0.001). Consequently, the difference in total body mass among the three feeding methods was not significant regardless of the feeding area.
However, the difference appeared in the variance; therefore, the feeding area was affected not by the total mass but by the variability in final growth.

Individual differences in growth
The individual difference in the body mass differed greatly, especially for Feeding C, the narrow-area feeding method (Fig 10). Fig 11 presents the body mass growth and feeding history of the individuals with the maximum and minimum body masses in Feeding C. The difference in body mass between the maximum and minimum individuals was small during the early days of the experiment (days 0-20; Fig 11). However, the difference between the two individuals increased gradually. As shown in Fig 11(b) and 11(c), the individual with the maximum body mass ingested feed well; in contrast, the amount of feed ingested by the minimum body mass individual was relatively small. This tendency continued until the final day.

Discussion
In this study, we developed a simulation model including a feeding behavior and a growth simulation. From the simulation results, the growth variability of the Feeding C group, which was the narrow-area feeding group, was larger than Feeding A, which was the widearea feeding group (Fig 8). This tendency was also reported by Jørgensen et al. [7], who conducted a rearing experiment involving standing and flowing water conditions. The feed clumped in the standing water tank and dispersed in the flowing water tank. As a result, the variability of growth in the standing water (clumped feed) was larger than that in the flowing water (dispersed feed). We demonstrated our hypothesis that the growth difference from feeding the spatial distribution was determined by individual fish behavior under the simulated conditions. The novelty of this study is that the growth of the fish can be explained by individual fish behavior.
The mean body mass in the simulation was compared with a past experimental result, and the body mass growth values agreed well with Morales et al. [19]. Similar variability in growth was confirmed in a past study [7]. Thus, we considered that the present simulation results were reliable. Growth history of the body mass of all 100 individuals using the three feeding methods, without performing randomization every day. The mean and standard deviation of the body mass and previous experimental results [19] are also indicated.
https://doi.org/10.1371/journal.pone.0280017.g008 A significant difference was observed between the variance of the body mass of each distribution on day 90 (Fig 10). Additionally, the skewness of narrow-area feeding was larger than that of wide-area feeding. These results suggest that the feeding distribution had an impact on the variability in growth, but not the mean body mass. We visualized the body mass growth and feeding history of the maximum and minimum individuals in Feeding C group; as shown in Fig 11, the individual difference gradually increased along with rearing day. The individual with maximum body mass ingested feed well. In contrast, the individual with minimum body mass did not feed well. Maximum swimming speed in the simulation was proportional to the total length of each individual; therefore, the swimming speed of a large individual was fast. These large individuals accessed and ingested feed more quickly than small individuals, and consequently the larger individuals grew more rapidly. The width of the feeding area was narrow in the Feeding C group, so the few individuals that arrived first to the feeding area Growth history of the body mass of all 100 individuals using the three feeding methods, with randomization performed every day. The mean and standard deviation of the body mass and previous experimental results [19] are also indicated.
https://doi.org/10.1371/journal.pone.0280017.g009 dominated feeding activity. This tendency was confirmed by the results shown in Figs 6 and 7, as feeding by the Feeding C group was completed faster than the Feeding A group. As a result, the individual difference in the Feeding C group was larger than the other feeding methods. The process of unequal growth was simulated, which will be useful for aquaculture management.
We now discuss the suitability of the simulation for aquaculture management. Feeding cost is critical when managing an aquaculture farm, so many studies have focused on the optimum feeding amount or methods for fish growth. In the past, an actual rearing test was the only way to determine the optimum feeding amount or method. However, this method is costly and time-consuming. In contrast, the proposed simulation confirmed growth before actual rearing. Therefore, the optimum feeding method and amount can be determined easily on a PC. Another advantage of the proposed simulation method is that detailed information on the reared individuals can be easily obtained and visualized. Exceptionally large or small individuals always appear during rearing. However, the reasons for this were not clearly understood because it is difficult to obtain an individual growth history. Therefore, we visualized the maximum and minimum body mass growth of individuals and their feeding history, and explored the reasons for growth in detail (Fig 11). Exceptionally large or small individuals have a negative effect on management of an aquaculture farm, because selling price is unstable. A fish farmer can more easily develop a management plan if individual growth history can be predicted. Exceptionally large or small individuals should be excluded for stable management, and predicting the timing is useful for farmers. Therefore, the proposed simulation will be useful for preliminary planning and future management without the need for an empirical method, and will aid to decision making by aquaculture farmers.
One limitation of this study was that some parameters, such as w n , the field of view, and excess feeding ratio, were determined by trial and error. For an accurate simulation, these parameters should be based on experimental results.
The behavioral parameters, w n , and field of view affected the schooling behavior directly. For example, a school will disperse when the value of w 1 (weight of "Separation") is large and gather when w 1 is small. A recent study proposed a method to measure 3D fish behavior using multi-stereo imaging [25]. In the future, 3D movement paths, including behavior according to feed presence, should be measured to understand interference by each individual of the target species. We believe that our qualitative results are not dependent on these parameters. However, these parameters should be determined from 3D movement paths in laboratory experiments for more accurate simulations.
The growth simulation method has another future task in this study. Here, growth was determined by the FCE [19]. Therefore, we estimated the body mass of rainbow trout accurately. However, this method can only be applied to rainbow trout in the same condition as those in the previous experiment. Therefore, the applicability of this simulation is limited. A growth simulation study using the DEB model has been reported [10]. The DEB model was proposed by Kooijman [8] and is based on metabolic theory, including energy and mass budgets for all living organisms at the individual level. The DEB model or another biological model should be used to simulate the rearing conditions without previous experimental rearing information. This should improve the versatility of the simulation.

Conclusions
In this study, we demonstrated that growth differences according to feeding spatial distribution were determined by individual fish behavior under the simulated conditions. Feeding cost is a critical problem when managing an aquaculture farm. Therefore, many rearing studies have been performed to understand the optimum feeding amount. However, conventional rearing tests are costly and time-consuming. As a result, operation of the actual aquaculture farm relies on the farmer's experience. Here, a rearing simulation test was easily conducted within a short time using the proposed simulation method, and different methods were attempted on a PC. The results demonstrated that optimal aquaculture operation can be determined without using an empirical method. This simulation is applicable only to rainbow trout under specific conditions. In the future, with the improvements mentioned in the Discussion, we believe that the proposed simulation method can be used as a decision-making tool in aquaculture and will contribute to efficient management of aquaculture farms.
Supporting information S1 File. Simulated growth and feeding history of Feeding A, without performing randomization every day.